Rayshader Exploration

output: html_document: toc: true toc_float: true df_print: paged code_download: true —

Introduction:

What is rayshader? Why is this so cool and fun? During this presentation, we will talk about the way to use rayshader, why we want to use rayshader and also how to change from 2D graphics to 3D graphics using rayshader.

How to install rayshader:

First, in order to use rayshader, we need to install some library

# install.packages("devtools")
devtools::install_github("tylermorganwall/rayshader")

## Skipping install of 'rayshader' from a github remote, the SHA1 (0dc8faa7) has not changed since last install.
##   Use `force = TRUE` to force installation

After we finished downloading the rayshader, we can try to see whether the library works or not.

library(rayshader)

If the library does not run, you shouold consider downloading XQuarz: https://www.xquartz.org/. Once you finish downloading XQuarz, hopefully, you can use Rayshader in your computer!

Why do we want to use rayshader?

• Supporting 2D and 3D data visualisation

• Directly transform the ggplot2 objects into 3D plot

• Generate beautiful topographic 2D and 3D maps

• Create a cinematic depth of field post-processing effect to direct the user’s focus to important regions in the figure

With all of these features, it will be wonderful to consider rayshader to create beautiful visualization for your graphics!

Diffferent available features for Rayshader:

Rayshader has a single function to generate 3D plots using ggplot2 objects:

plot_gg Takes a ggplot2 object (or a list of two ggplot2 objects) and uses the fill or color aesthetic to transform the plot into a 3D surface. You can pass any of the arguments used to specify the camera and the background/shadow colors in plot_3d(), and manipulate the displayed 3D plot using render_camera() and render_depth().

library(ggplot2)
ggdiamonds = ggplot(diamonds) +
  stat_density_2d(aes(x = x, y = depth, fill = stat(nlevel)),
                  geom = "polygon", n = 100, bins = 10, contour = TRUE) +
  facet_wrap(clarity~.) +
  scale_fill_viridis_c(option = "A")

par(mfrow = c(1, 2))

plot_gg(ggdiamonds, width = 5, height = 5, raytrace = FALSE, preview = TRUE)
plot_gg(ggdiamonds, width = 5, height = 5, multicore = TRUE, scale = 250,
        zoom = 0.7, theta = 10, phi = 30, windowsize = c(800, 800))
Sys.sleep(0.2)
render_snapshot(clear = TRUE)

Rayshader also detects when the user passes the color aesthetic, and maps those values to 3D. If both color and fill are passed, however, rayshader will default to fill.

mtplot = ggplot(mtcars) +
  geom_point(aes(x = mpg, y = disp, color = cyl)) +
  scale_color_continuous(limits = c(0, 8))

par(mfrow = c(1, 2))
plot_gg(mtplot, width = 3.5, raytrace = FALSE, preview = TRUE)

plot_gg(mtplot, width = 3.5, multicore = TRUE, windowsize = c(800, 800),
        zoom = 0.85, phi = 35, theta = 30, sunangle = 225, soliddepth = -100)
Sys.sleep(0.2)
render_snapshot(clear = TRUE)

Rayshader has seven functions related to mapping:

ray_shade uses user specified light directions to calculate a global shadow map for an elevation matrix. By default, this also scales the light intensity at each point by the dot product of the mean ray direction and the surface normal (also implemented in function lamb_shade, this can be turned off by setting lambert=FALSE.

sphere_shade maps an RGB texture to a hillshade by spherical mapping. A texture can be generated with the create_texture function, or loaded from an image. sphere_shade also includes 7 built-in palettes: “imhof1”, “imhof2”, “imhof3”, imhof4“,”desert“,”bw“,”unicorn".

create_texture programmatically creates texture maps given five colors: a highlight, a shadow, a left fill light, a right fill light, and a center color for flat areas. The user can also optionally specify the colors at the corners, but create_texture will interpolate those if they aren’t given.

ambient_shade creates an ambient occlusion shadow layer, darkening areas that have less scattered light from the atmosphere. This results in valleys being darker than flat areas and ridges.

lamb_shade uses a single user specified light direction to calculate a local shadow map based on the dot product between the surface normal and the light direction for an elevation matrix. add_shadow takes two of the shadow maps above and combines them, scaling the second one (or, if the second is an RGB array, the matrix) as specified by the user.

add_overlay takes a 3 or 4-layer RGB/RGBA array and overlays it on the current map. If the map includes transparency, this is taken into account when overlaying the image. Otherwise, the user can specify a single color that will be marked as completely transparent, or set the full overlay as partly transparent.

Rayshader also has three functions to detect and add water to maps:

detect_water uses a flood-fill algorithm to detect bodies of water of a user-specified minimum area.

add_water uses the output of detect_water to add a water color to the map. The user can input their own color, or pass the name of one of the pre-defined palettes from sphere_shade to get a matching hue.

render_water adds a 3D tranparent water layer to 3D maps, after the rgl device has already been created. This can either add to a map that does not already have a water layer, or replace an existing water layer on the map.

Also included are two functions to add additional effects and information to your 3D visualizations:

render_highquality renders in the scene with a built-in pathtracer, powered by the rayrender package. Use this for high-quality maps with realistic light transport.

render_depth generates a depth of field effect for the 3D map. The user can specify the focal distance, focal length, and f-stop of the camera, as well as aperture shape and bokeh intensity. This either plots the image to the local device, or saves it to a file if given a filename.

render_label adds a text label to the x and y coordinate of the map at a specified altitude z (in units of the matrix). The altitude can either be specified relative to the elevation at that point (the default), or absolutely.

And four functions to display and save your visualizations:

plot_map Plots the current map. Accepts either a matrix or an array.

write_png Writes the current map to disk with a user-specified filename.

plot_3d Creates a 3D map, given a texture and an elevation matrix. You can customize the appearance of the map, as well as add a user-defined water level.

render_snapshot Saves an image of the current 3D view to disk (if given a filename), or plots the 3D view to the current device (useful for including images in R Markdown files). render_movie Creates and saves a mp4 file of the camera rotating around the 3D scene by either using a built-in orbit or by using one provided by the user.

Plotting map and sphere_shade

library(rayshader)

#Here, I load a map with the raster package.
loadzip = tempfile()
download.file("https://tylermw.com/data/dem_01.tif.zip", loadzip)
localtif = raster::raster(unzip(loadzip, "dem_01.tif"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum unknown in Proj4 definition

unlink(loadzip)

#And convert it to a matrix:
elmat = raster_to_matrix(localtif)

#We use another one of rayshader's built-in textures:
elmat %>%
  sphere_shade(texture = "desert") %>%
  plot_map()

#sphere_shade can shift the sun direction:
elmat %>%
  sphere_shade(sunangle = 45, texture = "desert") %>%
  plot_map()

Detecting and Adding Water

#detect_water and add_water adds a water layer to the map:
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  plot_map()

Adding Shadow

#And we can add a raytraced layer from that sun direction as well:
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat), 0.5) %>%
  plot_map()

Adding Ambient Occlusion Shadow Layer

#And here we add an ambient occlusion shadow layer, which models 
#lighting from atmospheric scattering:

elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat), 0.5) %>%
  add_shadow(ambient_shade(elmat), 0) %>%
  plot_map()

Plotting in 3D

elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat, zscale = 3), 0.5) %>%
  add_shadow(ambient_shade(elmat), 0) %>%
  plot_3d(elmat, zscale = 10, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800))
Sys.sleep(0.2)
render_snapshot()

Adding Scale

render_camera(fov = 0, theta = 60, zoom = 0.75, phi = 45)
render_scalebar(limits=c(0, 5, 10),label_unit = "km",position = "W", y=50,
                scale_length = c(0.33,1))
render_compass(position = "E")
render_snapshot(clear=TRUE)

The problem set:

Below is our problem set about the Hong Kong housing price. We have given the datasets, the 2d Plots and also some information below.

library(readxl)
library(sp)
library(ggplot2)
library(png)
library(grid)
library(rayshader)
library(rgl)

## Here is the 3 datasets
district_name = read_excel('District Name.xlsx',1)
population = read_excel('Population Data.xlsx',1)
colnames(population) = c("Chi_name", "Population")

hkmap = readRDS("HKG_adm1.rds")

# The [reprocessing part
map_data = data.frame(id=hkmap$ID_1, Code=hkmap$HASC_1, Eng_name=hkmap$NAME_1)
map_data$Code = gsub('HK.', '', as.character(map_data$Code))
map_data = merge(map_data, district_name, by = 'Eng_name')
hkmapdf = fortify(hkmap)

## Regions defined for each Polygons

map_data = merge(hkmapdf, map_data, by="id")
map_data = merge(map_data, population, by = "Chi_name")
map_data$Population = as.numeric(map_data$Population)

# Here is the code to create a 2D map in Hongkong.
# Map
map_bg = ggplot(map_data, aes(long, lat, group=group, fill = Population)) +
  geom_polygon() + # Shape
  scale_fill_gradient(limits=range(map_data$Population), 
                      low="#FFF3B0", high="#E09F3E") + # Population Density Color
  layer(geom="path", stat="identity", position="identity", 
       mapping=aes(x=long, y=lat, group=group, 
                   color=I('#FFFFFF'))) # Boarder Color
  
map_bg = map_bg + theme(legend.position = "none", 
                        axis.line=element_blank(), 
                        axis.text.x=element_blank(), axis.title.x=element_blank(),
                        axis.text.y=element_blank(), axis.title.y=element_blank(),
                        axis.ticks=element_blank(), 
                        panel.background = element_blank()) # Clean Everything
map_bg

# Save as PNG
xlim = ggplot_build(map_bg)$layout$panel_scales_x[[1]]$range$range
ylim = ggplot_build(map_bg)$layout$panel_scales_y[[1]]$range$range
ggsave('map_bg.png', width = diff(xlim)*40, height = diff(ylim)*40, units = "cm")

# Real Estate Dataset
estate_df = readr::read_csv('https://raw.githubusercontent.com/cydalytics/HK_Properties_Price_Distribution/master/real_estate_master_df.csv')

## Rows: 195 Columns: 6

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): Real_Estate_Chi, Price_Per_SqFeet_Apr2020, Price_Per_SqFeet_Mar2020...
## dbl (2): Latitude, Longitude

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

estate_df$apr_price = as.numeric(gsub('[^0-9]', '', estate_df$Price_Per_SqFeet_Apr2020))
estate_df$mar_price = as.numeric(gsub('[^0-9]', '', estate_df$Price_Per_SqFeet_Mar2020))

# Read Background Image
hk_map_bg = readPNG('map_bg.png')

Here is the 2D map that we need to turn into 3D.

# 2D Plot
library(ggplot2)
library(grid)
estate_price = ggplot(estate_df) + 
  annotation_custom(rasterGrob(hk_map_bg, width=unit(1,"npc"), height=unit(1,"npc")), 
                    -Inf, Inf, -Inf, Inf) + # Background
  xlim(xlim[1],xlim[2]) + # x-axis Mapping
  ylim(ylim[1],ylim[2]) + # y-axis Mapping
  geom_point(aes(x=Longitude, y=Latitude, color=apr_price), size=2) + # Points
  scale_colour_gradient(name = 'Price per square foot (real)\n(HKD)', 
                        limits=range(estate_df$apr_price), 
                        low="#FCB9B2", high="#B23A48") + # Price Density Color
  theme(axis.line=element_blank(), 
        axis.text.x=element_blank(), axis.title.x=element_blank(),
        axis.text.y=element_blank(), axis.title.y=element_blank(),
        axis.ticks=element_blank(), 
        panel.background = element_blank()) # Clean Everything
estate_price

ggsave('estate_price.png', width = diff(xlim)*40, height = diff(ylim)*40, units = "cm")

After seeing the graph, the two questions we have will be:

1. How is the population distribution in Hong Kong? Is the estate price higher in densely populated district?

2. How to transfer a 2D plot to a 3D plot? Do you think creating a 3D plot will be a good idea in here? Why is that?